Single-Particle Dynamics in the Vicinity of the Mott-Hubbard Metal-to-Insulator 
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The single-particle dynamics close to a metal-to-insulator transition induced by strong repulsive 
interaction between the electrons is investigated. The system is described by a half-filled Hubbard 
model which is treated by dynamic mean-field theory evaluated by high-resolution dynamic density- 
matrix renormalization. We provide theoretical spectra with momentum resolution which facilitate 
the comparison to photoelectron spectroscopy. 
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I. INTRODUCTION 

Strongly correlated systems are notoriously difficult to 
understand. In particular, this is true when parameter 
regions are considered where the system changes its be- 
havior qualitatively, for instance where new excitations 
emerge. This is precisely the situation around the metal- 
to-insulator transition of electrons which repel each other 
strongly on a tight-binding lattice due to the Coulomb 
interaction. 

We consider a half-filled lattice with one electron per 
site on average. For strong on-site repulsion, it is ener- 
getically forbidden that two electrons occupy the same 
site. Because there is one electron per site on average 
this implies that there is precisely one electron per site. 
No electron can move so that the system is insulating. 
The remaining degree of freedom is the orientation of 
the electron spins. This freedom gives rise to collective 
magnetic excitations, for instance spin waves. 

For weak interaction, on the other hand, the system 
is a metal in the sense of a Fermi liquid. This is true 
if we neglect ordering phenomena, i.e., we focus on the 
paramagnetic phase. Hence this system is dominated by 
the quasiparticle excitations. No collective modes play a 
major role. 

In the vicinity of the mctal-to-insulator transition the 
Fermi liquid must change such that precursors of the col- 
lective magnetic excitations of the insulator emerge. Vice 
versa, the collective excitations in the insulator must ac- 
quire finite life times so that they broaden considerably 
in the corresponding spectra if we pass in the insulating 
regime towards the metallic one. Such behavior can be 
measured by inelastic neutron scattering for example. 

Hence it is of particular interest to understand the dy- 
namics of the important excitations in the vicinity of an 
interaction induced metal-to-insulator transition. In the 
present work, we focus on the single-particle excitations, 
that means we investigate the single particle propaga- 
tion. The investigation is done in the framework of the 



dynamic mean-field theory (DMFT) which maps the lat- 
tice model to an effective single site problem with a self- 
consistency condition.^'-' *' The single site is coupled to 
a bath of non-interacting fermions which can be repre- 
sented as a semi-infinite chain. 

The semi-infinite chain and its head where the interac- 
tion takes place can be tackled numerically by algorithms 
for one-dimensional systems. We will use density-matrix 
renormalization for this purpose. Since the dynamic 
properties have to be determined we use the dynamic 
density-matrix renormalization with correction vector. ' '^^ 
Thereby, we dispose of a numerical means which can re- 
solve spectral properties not only at low energies but also 
at high energies. ' This makes it possible to determine 
sharp spectral features away from the Fermi energy also 
quantitatively. ^" 

The main aim of the present article is to provide com- 
prehensive spectral data in the above mentioned theo- 
retical framework. In particular, we present momentum- 
resolved single-particle spectral densities in the insulating 
and in the metallic regime. 

In the following Sec. I A, we will present the model 
under investigation, i. e. the half-filled one-band Hub- 
bard model with semi-elliptic density of states (Bethe 
lattice) at zero temperature. Next, in Sec. II, we will dis- 
cuss the method. We use the dynamic mean-field theory 
(DMFT, Sec. II A) which maps the problem to an effec- 
tive single-impurity Anderson model (SIAM, Sec. II B) 
with a self-consistency condition. The SIAM is treated 
numerically by a dynamic density-matrix renormaliza- 
tion (D-DMRG) impurity solver (Sec. II C). The re- 
sults of this approach are presented in Sec. HI, with 
the spectral properties in the paramagnetic insulating 
phase (Sec. HI A) and the paramagnetic metallic phase 
(Sec. IIIB). The results reveal insights in the nature of 
the Mott-Hubbard metal-to-insulator transition and pro- 
vide a numerical proof for previous investigations which 
were based on the hypothesis of the separation of energy 
scales. Finally, a summary will be given in Sec. IV. 
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A. Model 

We use the generic model for strongly correlated sys- 
tems in solids, namely the Hubbard model. It is 
given by the Hamiltonian 

U = uY,{n,^^ - 1/2) (n,,i - 1/2) - t cl^c^^^ (1) 

where i, j denote sites on a lattice with (i, j) being near- 
est neighbors, a G {f , i} the spin and c^.^^ the electron 
annihilation (creation) operator and riijo- their occupa- 
tion operator. This model contains the basic ingredients 
of the problem, i.e., a local repulsive interaction and a ki- 
netic energy consisting of hopping from site to site. The 
interaction is diagonal in real space and hence tends to 
make the cigcnstates local in real space. The kinetic en- 
ergy is diagonal in momentum space and hence tends to 
make the eigenstates local in momentum space implying 
extended states in real space. So the interaction favors 
an insulating phase whereas the kinetic energy favors a 
conducting, metallic phase. They are the antagonists 
and depending on their relative strength the system is 
metallic or insulating. 

The transition between these two phases has been a 
long-standing issue. By now, evidence emerges that 
the transition is continuous at zero temperature^"'^'''^''' 
whereas it is of first order at any finite temperature.^'' 
The continuous transition at zero temperature is pecu- 
liar. It cannot be seen as an ordinary second-order tran- 
sition. It should rather be seen as a marginal first-order 
transition with zero hysteresis. At zero temperature, the 
free energy is continuously differentiable at the transition 
because the ground-state energy is continuously differen- 
tiable. The residual entropy jumps even at zero tem- 
perature because the ground-state of the metal is the 
non-degenerate Fermi sea. The insulator, however, is 
governed by the spin degrees of freedom. Since we do 
not consider any long-range order they remain free down 
to zero energy. The freedom to choose the spin locally t 
or I in the paramagnetic insulating phase implies a finite 
residual entropy of In 2 per site. But at zero temperature, 
the discontinuous entropy does not imply jumps in the 
other thermodynamic quantities. 



II. METHOD 

In this section we present the details of the methods 
that we used to tackle the model and to compute the 
desired single-particle dynamics. The key points are the 
dynamic mean-field theory, the single-impurity Anderson 
model and the dynamic density-matrix renormalization. 



A. Dynamic mean-field theory 

The basic idea of the dynamic mean-field theory is to 
consider the lattice problem under study in the limit of 
infinite coordination number z ^ 00. This means that 
the lattice is generalized in a way that z is a tunable 
parameter. Very often, this is realized by considering 
related lattices in d-dimensional space. Then z (x d and 
the desired limit is the limit of infinite dimensions. 

The limit z — > 00 is well-defined only for a particular 
scaling of the matrix elements in the Hamiltonian (1) 
which link different sites. In particular, the hopping to 
the nearest-neighbors must be scaled like 1/y/z.^'- Then 
the single-particle propagator between site j and site 
i scales like 

cx zll*-JII/2 (2) 

where we use — j|| for the taxi cab metric. This metric 
counts how many hopping processes are at least necessary 
to go from 7 to i. 

Considering the standard diagrammatic expansion in 
powers of the interaction one realizes that the scaling (2) 
implies that many diagrams vanish for z —^ 00. Only 
those which are either completely local or which are suf- 
ficiently numerous lead to non-vanishing contributions. 
Closer inspection shows that only those dressed skeleton 
diagrams lead to non-vanishing contributions which are 
completely local J-- This is called the collapse of the dia- 
grams because only a single site occurs. For the proper 
self-energy E this implies that it is local; only 'Su needs 
to be considered. 

In spite of the very substantial reduction of the number 
of relevant diagrams their summation is not directly pos- 
sible. Hence the dynamic mean-field theory requires an 
additional element. This is the observation that the same 
local dressed skeleton diagrams occur in the treatment of 
a single-impurity Anderson model (SIAM) where the in- 
teraction takes place at the impurity which is coupled 
to a fermionic bath. Hence, the same self-energy is 
obtained if the local dressed propagators are the same. 
Let us use capital letters for the local dressed propaga- 
tor G := Gii and for the local self-energy S := Si,; in 
the lattice model in DMFT. We use small letters for the 
local dressed propagator g and the self-energy a at the 
impurity of the effective single-site problem. Then the 
above statement reads 

J:{uj) = a{uj) (3a) 

if 

G{uj) = g{uj) (3b) 

holds. Note that this does not imply that G° = g° where 
we use the superscript ° for the bare, undressed propa- 
gators. 

For completeness, we mention that the self-consistency 
equations (3), which we derived from the expansion in 
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dressed skeleton diagrams, can also be obtained from the 
action by the so-called cavity argument. A derivation 
in the strong coupling limit is given in the review 18. An 
elegant derivation of the DMFT including the possibility 
to generalize from single site to cluster systems has been 
found by Potthoff'^'' in the framework of variational self- 
energy functionals. 

Eqs. (3) face us with the problem that wc do not know 
the bare propagator of the SIAM. This amounts up to 
a generic self-consistency problem where we first have to 
guess a and then we can verify whether the conditions 
(3) are fulfilled. 

Following this roadmap, two more fundamental rela- 
tions are needed, namely the Dyson equations of the lat- 
tice problem and of the SIAM. The local self-energy Yju 
acts on every site equally since wc consider the uniform, 
paramagnetic phase. For this reason wc could omit the 
subscript and pass on to S. This spatially constant en- 
ergy acts like a global energy shift so that it can be ac- 
counted for by 



is close enough to the previous assumption the iterations 
are stopped. Details of the precise meaning of how close 
is close enough are given in Sec. II C 4. 

The above sketched cycle can be used for any lattice, 
i.e., for any . For a particular with semi-elliptic 
spectral density (DOS) := -V'r3mG°(w) 



(7) 



G(uj) = G°{oj-Y.{uj)), 



(4) 



the cycle can be simplified considerably. A semi-elliptic 
spectral density p'^{u)) corresponds to the so-called Bethe 
lattice with infinite branching ratio. Note that we do 
not use any other feature from the Bethe lattice other 
than (7). Hence our calculation can be viewed as treat- 
ing a translationally invariant lattice with a semi-elliptic 
DOS. This is a good starting point for generic spectral 
densities since it is bounded and it possesses square-root 
singularities at the band edges like three-dimensional 
densities-of-states have. Otherwise it is featureless. 

The key feature of G° with semi-elliptic is its par- 

ticularly simple continued fraction representation with 
constant coefficients 



which is the Dyson equation for the lattice problem. 

In the single-site problem, the self-energy <j{ijj) acts 
only on the impurity site. The bare propagator from the 
impurity site to the impurity site is g^{'jj) so that the 
Dyson equation is given by the geometric series 



^0 



.9°(^) 



0\ri 



l^g^{u)a{ujy 
This relation can also be written as 



(5) 



(6) 



The above equations are used to set up an iteration 
cycle to find a self-consistent solution. It is illustrated 
in Fig. 1. Starting in the upper left corner, an initial 
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a = Z 



Figure 1: Schematic representation of the iterative self- 
consistency cycle for an arbitrary lattice. In the step marked 
with * the dynamics of an effective single-impurity Anderson 
model is calculated by a D-DMRG impurity solver. 

guess is used to define the SIAM which is then solved by 
an impurity solver providing a. Then a is used as lattice 
self-energy S providing via the lattice Dyson equation (4) 
the dressed propagator for the lattice. By the condition 
(3b) the dressed propagator of the SIAM is known which 
in turn yields via (6) its bare counterpart so that wc arc at 
the beginning of the cycle again. If the bare propagator 



G°(c.) 



uj - 



74 



(8) 



UJ ■ 



In order to use this fact for simplification wc write the 
bare propagator of the SIAM with the help of the so- 
called hybridization function r(cL') as 



1 



cij — r(cij) 



(9) 



where the continued fraction of r(w) shall be 
parametrized like 



y2 

F(.)^^ ^ 



(10) 



7? 



Then, based on Eqs. (6,3a) the dressed propagator of the 
SIAM reads 



1 



3(w) 



= cij — r(cij) — s(cj). 



(11) 



By Eqs. (4,8) the dressed propagator of the lattice can 
be expressed to be 



1 



G{u) 



UJ — S(cij) 



^74 



w — s 



D 



74 



— E 



'74 



w — E ■ 



- S(u;) - (^74) G(c^), 



(12) 



where we exploited that the continued fraction docs not 
change when it is evaluated at a deeper level because its 
coefficients are constant. 
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Based on the self-consistency condition (3b) we set 
Eqs. (11,12) equal and obtain the simpler self-consistency 
condition 



(13) 



This equation is more easily used because it provides a 
direct way to compute the hybridization function T(uj) of 
the next iteration of the SIAM from the dressed lattice 
propagator G. No explicit computation of intermediate 
self-energies is needed which can be severely hampered 
by numerical errors, see below. A first conclusion which 
can be drawn from (13) is that V = D/2. 

The resulting cycle is given graphically in Fig. 2. This 
form of the iteration cycle is used in the present work. 



g 



initial 
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self-consistency? 



g 



(13) 



g 

(3b) 
G = g 



Figure 2: Schematic representation of the iterative self- 
consistency cycle for semi-elliptic density of states (DOS) 
with nearest-neighbor hopping. The step * stands for the 
D-DMRG impurity solver. 



B. Single-Impurity Anderson Model 

There are various ways to set up the effective single- 
impurity problem because the precise topology of the 
fermionic bath, to which the impurity couples, does not 
matter. As we have seen in the equations above we are 
only concerned with the local quantities on the impurity. 
Hence the only relevant property of the fermionic bath is 
the bare propagator or the hybridization function F, 
respectively. 

Because we aim at the application of density-matrix 
renormalization, which has been developed for one- 
dimensional systems in particular, we favor the represen- 
tation of the bath as semi-infinite chain, see Fig. 3. This 
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Figure 3: (Color online) Sketch of the semi-infinite chain 
which is used as topology for the single-impurity Anderson 
model. 

is no restriction since mathematically any hybridization 
function which has a non-negative spectral density can 
be represented as continued fraction"" as given in Eq. 
(10). By complete induction it is straightforward to see 
that the choice of the hopping elements shown in Fig. 3 
implies Eq. (10). 



The Hamiltonian of the SIAM represented as linear 
chain reads"^^ 



n = U (nd^ - 1/2) (nd^t - V2) + VY, (4co,<, + h-C 

(T 

OC 

+ 7n(4.c„+i,^+h.c.). (14) 



n=0, (7 

The correlations stem from the density-density coupling 
on the impurity site {S^^ operators) at the head of a semi- 
infinite chain of free fermions (c^^^ operators). For com- 
prehensive information on the physics of single-impurity 
Anderson models the reader is referred to Ref. 21 . 

The determination of the local propagator g{ijj) in the 
effective SIAM (14) for a given hybridization function 
Tiuj) is the most difficult part in the self-consistency cy- 
cles (cf. Figs. 1 and 2). The necessary tool is called the 
impurity solver. 

The dynamics one has to compute is the dynamics of 
the fermionic single-particle propagator of the impurity 
electron. Aiming at the properties at T = it reads 



Wo 



1 



di 



uj + bj - {H - En)'^'' 



(15) 



Here the ground-state is denoted by |0) and its energy by 
Eq. Since we aim at the paramagnetic phase the prop- 
agator has no dependence on the spin index a and we 
do not denote it. The frequencies uj and rj are real. The 
standard retarded propagator is retrieved for 77 —>■ 0+ 



gniui) = lim giiv + iij). 



(16) 



The wanted quantity is the spectral density p{uj) := 
—n^^3mgYi{uj). If necessary the real part of gniuj) can 
be obtained from p{uj) by the Kramers-Kronig relation 



00 

megniuj)^V J 



(17) 



where V denotes the principal value. 



C. Dynamic Density-Matrix Renormalization 

1. Algorithm 

Because there are essentially no analytic approaches 
available which solve the general SIAM, reliable numeri- 
cal approaches for dynamic correlations are very impor- 
tant. Our choice is the dynamic density-matrix renor- 
malization which is an excellent tool for calculations for 
open one-dimensional systems at zero temperature. 

To obtain g{uj) for a given F(a;) we use a combination 
of the dynamic density-matrix renormalization '■''^'^^ 
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(D-DMRG) and the least-bias deconvolution 
algorithm.-' The D-DMRG is the generalization of 
the standard density-matrix renormalization group^'*' 
(DMRG) for the calculation of dynamic quantities. A 
key advantage over other T = approaches like the 
numerical renormalization group is that the resolution 
can be high not only at low energies but also at high 
energies.'* We implemented the D-DMRG in a correction 
vector (CV) scheme' with optimized direct matrix 
inversion, for details see Refs. 9 and 25. 

The actual calculations are not performed for the 
fermionic representation (14) but for a spin represen- 
tation involving two semi-infinite spin chains which are 
coupled at their heads. This spin representation re- 
places each fermionic level with states empty and occu- 
pied by the two spin states 1 and |. The change of ba- 
sis between both representations is accomplished by two 
Jordan- Wigner transformations, one for the "f electrons 
and one for the | electrons, see Ref. !). The resolvents 
appearing in Eq. (15) are expressed as resolvents of the 
equivalent spin system involving spin flips at the head of 
the chains. The advantage for the DMRG approach is 
that in the spin chain description each site carries only 
a two dimensional Hilbert space instead of a four dimen- 
sional one. 

The correction vector D-DMRG provides data points 
at given frequencies lo ~ uji for finite values of r] = rji 



1^ . 
Jm5(cjj + 1?7, 

TT 

p{uj) (g) L^^ {uj] 

oo 

/■ riip{uj)du! 

TT J {uj - UJi)'^ + rjf 



(18a) 
(18b) 

(18c) 



where OL,,; (lu) stands for the convolution with a 
Lorentzian of width rji. No data can be obtained directly 
at ry = since the inversion of the Hamiltonian is singular 
and cannot be achieved numerically in a stable way. Fur- 
thermore, we cannot numerically treat an infinite chain 
but only a finite one. 

Henceforth, we will call data {gi} raw data because it is 
obtained at finite values of rji . Apart from the constraint 
from matrix inversion, a sufficiently large broadening rji 
is used in order to calculate continuous spectral densities 
resembling the spectral densities in the thermodynamic 
limit. The finite-size effects are washed out. In a sec- 
ond step, we aim at retrieving the unbroadened spectral 
density p{uj) as well as possible. 

We use the D-DMRG raw data calculated for a fi- 
nite set of frequencies uji with non-zero broadenings rji 
as input of a deconvolution scheme to extract the in- 
formation on the spectral density p{i-u), i.e. the rele- 
vant dynamic properties of the infinite system. The de- 
convolution scheme of our choice is the least-bias (LB) 
approach" ' which belongs to the class of maximum en- 
tropy methods.-'' The combination of D-DMRG and LB 



deconvolution has already proven itself in the investiga- 
tion of the high-energy dynamics of the SIAM by provid- 
ing a well-controlled resolution for all energies."''" ' By 
construction, the LB ansatz yields a positive and con- 
tinuous result for the density of states p{oj) (DOS). The 
positivcness assures the causality of the solution, which is 
necessary for calculating the continued fraction, see App. 
A, whereby the DMFT self-consistency cycle is closed. 
By using the LB deconvolution we avoid any kind of ar- 
bitrariness introduced by the choice of a discretization 
mesh. 

Of course, the direct computation of a continuous 
density from data obtained for a finite system involves 
an approximation. In practice, one must pay attention 
whether the system size and the broadening considered 
make it possible to deduce continuous densities, see also 
below. 

Nishimoto and co-workers"''"'^ advocate a different ap- 
proach which has been baptized "Fixed-Energy" DMFT 
where the problem of deducing densities is circumvented 
by discretizing the continuous frequency in a number of 
bins. For these bins the spectral weights are computed. 
On the one hand, spectral weights are better behaved for 
finite system sizes. But on the other hand, the DMFT 
self-consistency holds only for the continuous quantities 
in Eq. (3). So any extension to discrete quantities in- 
troduces some ambiguity.""' Balancing these two contrary 
aspects we have chosen to work with the continuous spec- 
tral densities as obtained from the LB ansatz. 

In the framework of the D-DMRG Jeckelmann has for- 
mulated a variational principle for the wanted dynamic 
correlation. '" A certain functional is minimum for the 
correct dynamic correlation. Its value yields the gi as 
defined in (18a). Due to the minimum principle this ap- 
proach is rather robust. Its main advantage is that even 
if the numerical dynamic correlation possesses an error d 
the values gi will be exact up to quadratic order . From 
the algorithmic view point the minimization required for 
Jeckclmann's approach is more demanding than the ma- 
trix inversion in the D-DMRG with correction vector. 
Hence we stick to an optimized matrix inversion.'''" ' In 
our opinion, the optimum strategy is to search for the 
best correction vector by optimized matrix inversion and 
then to use the variational functional for the evaluation 
of the gi. 



2. Deconvolution and Finite-Size Effects 

The choice of the rji for a given set of uji is restricted 
by the considered chain length L. If the r]i are chosen too 
small, the extracted spectral properties reveal too strong 
finite-size effects. If the r]i are chosen too large, essen- 
tial features in the spectral properties cannot be resolved 
properly and therefore remain smeared out, even after 
the deconvolution. Thus, the optimum value of each 77^ 
has to been chosen with care depending on the width of 
the features which have to be resolved. In our approach 
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the optimum choice is characterized by 

7?,; « Alu, > W*/L (19) 

where Auji denotes the distance between two neighboring 
frequencies and W* is the bandwidth. The approximate 
equahty between rji and Aw,; reflects the fact that it is not 
reasonable to numerically measure the spectral density 
which is broadened by 77.^ at frequencies much closer than 
T]i. This does not lead to an increase of relevant informa- 
tion about the underlying unbroadened spectral density. 
In contrast, it might happen that small errors lead to 
inconsistencies of the too closely positioned data points 
so that the subsequent deconvolution leads to strongly 
oscillatory results.-' 

The inequality Aui » W* / L in (19) results from the 
approximation of the infinite system under study by a 
finite system. If the points tOi, at which the broadened 
spectral density is measured numerically, were too closely 
spaced the finite-size effects would be discernible. Hence 
the raw data {gi} could not be processed as if it resulted 
from a continuous spectral density. 

Finally, we point out, that apart from the resolution 
aspect, the choice of rji influences the stability and the 
performance of the correction vector D-DMRG. For very 
small values of r]i the matrix inversions are numerically 
very demanding and time-consuming. 



3. Continuous Fraction 

The LB algorithm and the Kramers-Kronig transfor- 
mation allow us to extract the full unbroadened Green 
function g{uj) to good accuracy from the raw D-DMRG 
data for an arbitrary to. Thus we can have an arbi- 
trary number of supporting points for further numerical 
calculations. The hopping elements 7„ of the effective 
SIAM (14) for the next iteration in the self-consistency 
cycle are calculated recursively by a continued fraction 
expansion.-"' Note, that even tiny regions of small neg- 
ative spectral density spoil the continued fraction expan- 
sion completely. The details of the calculation of the 
continued fraction are given in App. A. 



4-. Numerical Criteria for Self-Consistency 

We have realized the iterative cycle shown in Fig. 2. 
Since we perform the calculation of the spectral den- 
sities numerically, including an advanced deconvolution 
scheme, the calculation of two identical spectral densities 
within the iterative self-consistency cycle is not possible 
by principle. Consequently, a certain error tolerance has 
to be included to the termination conditions. 

A solution within our DMFT self-consistency cycle is 
understood as self-consistent, if at least the following con- 
ditions are fulfilled: 



1. Two deconvolved spectral densities calculated con- 
secutively within the self-consistency cycle obey for 
all Lu, 

|p«(c.)-p('+i)(c.)|<e, (20) 

where e = 4 x 10^^ /D for the metallic and e = 
1Q~^/D for the insulating solutions. 

2. The average double occupancy and the ground- 
state energy per lattice site calculated in succession 
within the self-consistency cycle are stable. That 
means that they differ by less than 10^^%. 

For the insulating solutions, we additionally require that 
the static single-particle gaps calculated consecutively by 
means of the standard DMRG algorithm differ by less 
then 1%. The higher value of e for the metallic solutions 
has to be chosen because the metallic spectral densities 
display a high complexity including many sharp features. 
In contrast, the spectral densities of the insulating solu- 
tions are rather featureless and thus converge more easily. 

We stress that the upper conditions are the minimum 
requirement and that most of our solutions surpass them 
considerably. But we observed that it is sufhcient to ter- 
minate the self-consistency cycle once the upper condi- 
tions are fulfilled. The physical quantities obtained from 
the continued fractions computed in two consecutive it- 
erations are almost identical. Further iterations yield 
no further improvement. Moreover, we have confirmed 
that our final continued fraction coefficients do not dis- 
play persistent trends of change if iterated further. They 
rather display very small oscillations around the results 
which we consider to be converged. 

5. Numerical Stability of the Metallic and the Insulating 
Solution 

There is a parameter region in the Hubbard model in 
infinite dimensions with an interaction U between Ud 
and Uc2 where both solutions are locally stable. In a 
numerical treatment of this situation one has to introduce 
a bias, preferably minute, which selects the solution we 
want to find. 

This can be done elegantly by exploiting an even-odd 
effect. If we consider only the impurity site, the bare 
DOS p'^{uj) consists of a single S-peak at a; = 0. Any 
arbitrary finite repulsion U > splits this peak into two 
so that a gap A ~ U/2 occurs. The corresponding self- 
energy has a pole at w = which is necessary to induce 
the splitting of the central peak. Hence this self-energy 
is the self-energy of an insulator. 

Any SIAM with an odd number of fermionic sites (one 
impurity site and an even number of bath sites) has a bare 
DOS p°(w) with a (5-peak at a; = 0. This lowest energy 
level can be occupied or unoccupied without changing 
the ground-state energy. A finite repulsive interaction 
will lift this degeneracy partly and split the central S- 
peak in two and a gap occurs. Of course, if the levels in 
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the bare system lie very close to each other the occurring 
gap can be very small. But in any case, the solution 
resembles an insulating one. We conclude that the use of 
a SIAM with an odd number of sites implies a small bias 
towards an insulating solution. Obviously, the relative 
bias is of the order of the inverse system size l/L. 

Conversely, considering a SIAM with an even number 
of sites amounts up to a small bias towards a self-energy 
belonging to a metallic phase. The imaginary part of E 
at w = is zero. 

While the metallic solutions correspond to the Fermi 
sea as unique ground-state the paramagnetic insulating 
solution implies a macroscopically large degeneracy in the 
lattice problem. In the SIAM, the degeneracy is only two. 
It stems from the fact that the spin of the impurity is not 
completely screened by the coupling to the bath because 
the bath is gapped by A. So there is no DOS at low ener- 
gies and the renormalization of the unscreened moment 
stops at the lower cut-off A. Indeed, we found that the 
insulating baths show a two-fold degenerate ground-state 
or two almost degenerate low-lying states. 

The spin degeneracy in the insulator poses a problem 
to the numerical treatment. It happens easily that the 
DMRG algorithm looking for the ground-state provides 
a state with a certain finite magnetization in z-direction. 
If such a solution is used further in the iteration cycle 
results are produced which do not describe the paramag- 
netic insulating phase. In order to avoid this problem, we 
first search for the ground-state |1) with energy Ei using 
a standard Lanczos or Davidson algorithm starting from 
an initial guess The degenerate or almost degener- 
ate second state |2) is found from a second run starting 
from the orthogonal initial guess I/2) = — 
With both states we can construct a ground-state with- 
out magnetization. We use the normalized superposition 
|0) = ai|l) + a2|2) with 1 = |qiP + |a2p where the 
coefficients are chosen such that 

(0|d,4|0) = 1/2. (21) 

for <T = {t, i} which ensures the absence of a magnetiza- 
tion in z-direction. 

Finally, we note that the treatment of the spectral den- 
sity in the region of the single-particle gap A in the insu- 
lating phase requires care. The LB deconvolution scheme 
does not allow for zero spectral density so that some ar- 
tifacts need to be removed before the continued fraction 
expansion can be calculated as explained in Appendix A. 
The details of the removal of deconvolution artifacts are 
given in App. B. 

III. RESULTS 

Here we present the results for various spectral densi- 
ties which we have obtained by the approach described 
so far. First we focus on the insulating solutions; then 
we pass on to the metallic solutions. The results are in- 
terpreted as well as compared to previous results. 



A. Insulator 

1. Local Spectral Densities 

The computation of the insulating spectral densities 
is performed for 121 fcrmionic sites, including the impu- 
rity. Larger system sizes provide identical solutions (not 
shown). This is related to the relatively low complexity 
of the insulating solutions which do not display particular 
features besides the insulating gap. 

We calculate the raw D-DMRG data on a grid using 
mostly lS.LOi = 7;^ = O.ID. Close to the edges of the 
Hubbard bands we increase the resolution by choosing 
Aw,; = rii = 0.05D and partially 0.025D or even Awi = 
T]i = OMD for U = 2AD. This is done to resolve the 
gap and the band edges properly. Note, that such a small 
broadening is useful here, although rji = O.OID does not 
obey the inequality (19), because it helps to decide where 
the spectral density is finite or zero. 
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Figure 4: (Color online) Spectral densities p{u)) in the insulat- 
ing regime (solid line) calculated with a chain of 121 fermions. 
The dashed lines (almost coinciding) show the results from 
strong coupling perturbation theory. '~ 

The insulating solutions in Fig. 4 clearly show the lower 
and the upper Hubbard band each with a band width 
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of nearly 2D. The bands are separated by a well pro- 
nounced gap of 2 A:''' The center of each Hubbard band 
is located at roughly uj = ±'^/2. The shape of the density 
in the Hubbard bands in our insulating solutions is essen- 
tially smooth and reveals no significant features. They 
agree very well with the perturbative results by Gebhard 
et al. '" which are based on a strong coupling expansion 
in 1/(7. Only for very small gaps A w 0.2D a noticeable 
deviation is discernible, see Fig. 5. Clearly, these devia- 
tions result from the limited number of terms which arc 
available in the perturbation series. 
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Figure 5: (Color online) Data from Fig. 4 on a zoomed scale 
close to the inner band edges. 

In Ref. 33 a DMFT calculation based on another 
DMRG impurity solver, the Lanczos method, showed 
many substructures in the Hubbard bands. The com- 
parison to our featureless bands leads us to the conclu- 
sion that the peaky substructures in Ref. 33 stem from 
finite-size effects. First, such effects enter because only 
45 fermionic sites were considered. Second, the reliable 
application of the Lanczos impurity solver would require 
a large number (sa 100) of target states which cannot be 
handled by the DMRG. 

The results in Ref. 3)4 differ from ours only in one im- 
portant aspect. Nishimoto et al. found an upturn of the 
spectral density p{lS) close to the inner band edge for 
small gaps which is completely absent in our data. We 
believe that this difference stems from the different treat- 
ment of the spectral density in the gap region, for our 
approach see Appendix B. In our calculations, Nishi- 
moto's approach to use the extrapolated static gap as a 
cutoff for the spectral density did not lead to stable self- 
consistent solutions. On iteration, the spectral density 
becomes peakicr and peakier. The normalization does 
not hold anymore. 

But we found that a slightly too large gap value tends 
to produce features similar to the upturns in Ref. 34. 
Our approach is corroborated further by the excellent 
agreement between our value of the interaction C/ci = 
(2.38 ± 0.02)D below which no insulating solution exists 



and the one obtained by extrapolating the i/(7 expansion 
of the ground-state energy. The latter is supported also 
by quantum Monte Carlo results. In spite of the evi- 
dence mentioned above, one may certainly argue in favor 
of the Fixed-Energy approach used in Ref. 34 or in fa- 
vor of the approach used here because both numerical 
approaches comprise some approximations. For further 
comparison, see also the next subsection Sec. HI A 2. 

The perturbation theory'" implies a square-root be- 
havior at the band edges. The extremely good agree- 
ment between the perturbative and the numerical results 
in Figs. 4 and 5 suggests that a square- root onset is the 
correct description 

PuHB(c^)(x(t^-A)i/2 (22) 

at the inner band edges of the upper Hubbard band 
(UHB). 

In order to investigate the behavior of the spectral den- 
sity quantitatively near the band edges, see Fig. 5, we fit 

/(c^) = e(cc>-A)[(cc;-A)"P3(^-A)] (23) 

to the inner side of the upper Hubbard band (not shown). 
Here Q(x) denotes the Heaviside function, Pnix) is an 
nth order polynomial 

n 

P„(x)=^c,x^ (24) 

1=0 

and a, A, and the c,; are the fit parameters. The param- 
eter a represents the leading exponent of the power law 
onset. Note, that a perfect square-root onset cannot be 
extracted using the LB deconvolution scheme since the 
ansatz for the unbroadened data is built from exponen- 
tial functions.-' ' Hence the fit results depend on a number 
of assumptions, for instance the degree n of Pn{x) and 
the interval in which the least-square fit is done. Fig. 6 
summarizes our results; the ambiguities are accounted 
for by the error bars. Our results excellently confirm 
the square-root onset. The agreement is very good for 
large values of the interaction, i.e., deep in the insulating 
regime. Closer to the closure of the gap the uncertainties 
grow. 

The outer band edges of the Hubbard bands do not 
show any sharp decay. Instead, they reveal a smooth 
drop for — s- cxd which becomes even smoother as U 
decreases, similar to the behavior of the metallic solu- 
tions, cf. Sec. IIIB. 

2. Insulating Gap 

The gap A is the characteristic quantity of the insu- 
lator. Because of the impossibility to determine it di- 
rectly from the LB deconvolved spectral density we use 
the value which is found from the fit (23) with fixed ex- 
ponent a = 1/2. We stress that this determination of the 
gap A does not enter the iterations of the self-consistency 
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Figure 6: (Color online) Leading exponent a in the behavior 
of the spectral density near the band edges {uj —r A) for 
various values of U. The dashed line marks the value a = 0.5 
suggested by perturbation theory. 



cycle, except for the removal of artifacts, see App. B 2. It 
is done only at the end of the iterations once the solution 
can be considered to be self-consistent according to the 
criteria in Sec. II C 4. This is in contrast to the Fixed- 
Energy procedure chosen in Ref. 34 where the gap enters 
in the self-consistency explicitly. Again, it is a priori not 
clear which approach is more advantageous. All we found 
was that we could not reach stable self-consistency when 
using extrapolated gaps in the devonvolution. 
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perturbation theory expanding in powers of ^/u?'^ 

Closer to the vanishing gap the behavior of A([/) is 
not linear anymore but displays a significant curvature. 
Qualitatively, this is similar to what has been observed 
previously in iterated perturbation theory " (IPT), by 
the local moment approach "^ (LMA), or by the Lanczos 
approach D-DMRG." This observation leads us to fit a 
power law 

/a(C/) = ([/ - C/cO'^K + a^iU - C/ei)] (25) 

with the fit parameters Oj, Ud, and C, which denotes 
the leading exponent of A(C/) close to Uci- From this 
fit f/ci = (2.38 ± 0.02)D and C = 0.72 ± 0.05 is found. 
Thus, our value for Ud agrees excellently with most of 
the previous results, see Tab. I. 



3. Self-Energy 

The self-energy E is not required in the self-consistency 
cycle shown in Fig. 2, but it is accessible via Eqs. (3a, 6). 
By means of the self-consistently determined full propa- 
gator G{lo) and the bare propagator g^{Lo) of the SIAM 
we compute 



S(^) 



1 



1 



G{ujy 



(26) 



Note, that this way to access the self-energy can be prone 
to numerical problems. If the self-energy is a small quan- 
tity it can acquire a substantial relative error if it is 
computed according to (26) as difference of two large 
quantities. In the extreme case, it can happen that the 
imaginary part of the retarded self-energy computed via 
(26) is positive, i.e., it displays the wrong sign and vi- 
olates causality. This will be a problem in the metallic 
solutions at weak coupling, see Sec. IIIB. 

The indirectly calculated self-energy in the insulating 
regime displays the expected features, sec Fig. 8. Except 
for a i5-peak in the imaginary part at a; = 0, S(cl;) is 
strictly real inside the gap 



3mE(cj) = —cit5{uj) for |w| < A , 



(27) 



2.5 



3 

U/D 
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with c > 0. Correspondingly, the real part of S(toi) shows 
an l/w behavior: d\tTi{Lo) = c/ui + O{uo). 



Figure 7: (Color online) Single-particle gap A from the insu- 
lating solutions and the pseudogap from the metallic solutions 
vs. the interaction U . The determination of the pseudogap is 
explained in Sec. IIIB. For comparison the results from per- 
turbation theory (PT) up to 0{U^^) in the insulating regime 
are also shown." ' 

The results are plotted in Fig. 7 as a function of U . 
As expected, the gap A decreases continuously with de- 
creasing U until it vanishes at a critical value Ud- For 
U > 3.4D our A([/) agrees well with the results from 



4- Momentum-Resolved Spectral Densities 

So far we discussed the local spectral density p{lo) 
which is related to the imaginary part of the local prop- 
agator Gii{uj). This is the quantity which matters in 
the self-consistency of the DMFT approach. But once a 
self-consistent solution is found one can go a step further 
towards spatial correlations. While the self-energy and 
the skeleton diagrams are local in DMFT (see above) one- 
particle reducible quantities like the propagators are not. 
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w/D 

Figure 8: (Color online) Real and imaginary part of the self- 
energy 'S{uj) in the insulating regime. The 5-peak in the imag- 
inary part (lower panel) is not shown. 

Hence it makes sense to discuss Gij . Because we study a 
translationally invariant system we study the propagator 
in momentum space G'fe(ijj).*'" Its Dyson equation reads 



By varying the bare fermionic energy we can assess the 
effect of different momenta. In particular, if we change ek 
from negative to positive values we can view the change of 
the single-particle response on passing from the dynamics 
of a hole to the one of an electron. 

Figs. 9 and 10 display the spectral density of Gk{oj) 

p{Lj,ek) ■.= --JmGk{uj). (29) 

TT 

Tiny wiggles in the intensity as function of the frequency 
are numerical artifacts of the self-energy computed via 
(26). The gap between the lower and the upper Hub- 
bard band is clearly visible. In each Hubbard band there 
is a ridge of high density running from the upper right 
to the lower left edge. Neglecting the strong scattering 
we can interpret this ridge as the dispersion of a moving 
double occupancy in the upper Hubbard band. Corre- 
spondingly, there is a moving empty site in the lower 
Hubbard band. Since the band width of each Hubbard 
band is approximately equal to the band width of the 
bare dispersion we deduce that the hopping of the dou- 
ble occupancy or the empty site, respectively, equals the 
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Figure 9: (Color online) Spectral densities p{uj,ek) in the in- 
sulating regime for U = 2.4_D (top), U — 2.6D (middle), and 
U = 2.SD (bottom). 



bare hopping. This fact has been exploited also in the 
perturbative treatment.^'' 

Note, that the first approximation to this problem ob- 
tained by Hubbard^" yields only half the band width 
for each Hubbard band. But the famous improved solu- 
tion in Hubbard's third paper displays the correct band 
width for strong coupling. This treatment also includes 
scattering processes which take the strong scattering of 
the additional fermion into account. Hence it is very in- 
structive to compare our numerically exact result with 
the approximate solution in Ref. 39. Fig. 11 shows two 
examples for an almost vanishing gap at U = 1.8D and 
for a large gap a,t U — 4.0D. The similarity to the ex- 
act results in Fig. 9, upper panel, and in Fig. 10, lower 
panel, is striking. Both results, exact and approximate, 
show strongly scattered, hence broadened, excitations in 
the Hubbard bands. The ridges behave very similarly 
in shape and height. In particular for the large gap the 
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Figure 10: (Color online) Spectral densities p{u,ek) in the 
insulating regime for U = 3. 2D (top), U — 3.6D (middle), 
and U = 4,0 (bottom). 



results agree well. 



For the smaller value of the insulating gap. the most 
important difference is that the value of U differs. But 
if we accept that the interaction value in Hubbard's ap- 
proximation needs to be renormalized in a certain way 
the spectral densities behave qualitatively very similar. 
It must be stressed, however, that the results are quanti- 
tatively different. Comparing the upper left panel in Fig. 
4 with the corresponding result in Ref. 39 one notes that 
the spectral weight in the exact solution is found further 
away from the Fermi level at w = while in Hubbard's 
solution it is shifted towards a; = 0. 
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Figure 11: (Color online) Spectral densities p{u>,ek) in the 
insulating regime according to Hubbard's approximation in 
Ref. :'>!) for U = 1.8D (upper panel) and U = 4.0D (lower 
panel) 



B. Metal 



Local Spectral Densities 



Figure 12 shows the self-consistently calculated spec- 
tral densities p{u!) well in the metallic regime (0 < 
U < 2D). Figure 13 shows the corresponding densities 
p{uj) closer to the transition to the insulating solution at 
U = Uc2 ~ 3.0D. We used 160 fermionic sites (includ- 
ing one site representing the impurity) in the D-DMRG 
impurity solver. This implies a tiny bias towards the 
metallic solution, see Sec. II C 5. The raw data is ob- 
tained on an equidistant frequency grid with step size 
Aoji = O.ID using a constant broadening 77 = O.ID. It 
is deconvolved using the LB algorithm.- ' The spectral 
densities in Figs. 12 and 13 qualitatively display many 
of the features known from previous analytical^'*'^" and 
numericah'^''^'^'"'-'^'"'^'"''^'"'"''"''''"''''"^ ' investigations. 

Starting from the semi-elliptic DOS (7) at U ~ the 
central peak becomes narrower and narrower on increas- 
ing U. The diminishing width is the Kondo energy scale 
Tk. The weight in the central peak is the quasiparticle 
weight Z which we will discuss below. 

The height of the central peak does not change. The 
quasiparticle peak is pinned at the Fermi energy to its 
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Figure 12: (Color online) Spectral densities p{uj) in the metal- 
lic regime calculated with a chain of 160 fermionic sites. The 
dashed lines show results of a DMFT approach based on the 
numerical renormalization group (NRG) as impurity solver.""' 



non-interacting value 



p{uj = 0) = p"{uj = 0) 



2 



(30) 



as required by Luttinger's theorem for a momentum in- 
dependent self-energy.- As the DOS at zero frequency is 
not pre-determined in our approach, the pinning crite- 
rion (30) provides a sensitive check for the accuracy of 
the numerical solutions in the metallic regime. This cri- 
terion is very well satisfied, see Fig. 12 and the upper 
panels in Fig. 13. A maximum relative error of 3% oc- 
curs for interaction values close to the transition to the 
insulator, see lower panels in Fig. 13. The fulfillment of 
the pinning is an important prerequisite for an accurate 
determination of the quasiparticle weight Z below. 

Beyond U ^ D, the DOS shows side features. They 
are just shoulders which develop into the lower and upper 
Hubbard bands at larger interaction. As expected, they 
appear roughly at w = zLU/2. At J7 « 2D the Hubbard 
bands are well separated from the quasiparticle peak by 
a precursor of the gap A in the insulator: a pseudogap 
is formed. This pseudogap is determined from U = 2AD 



3 
Q 



0.5 



' 1 ' 1 


' 1 ' 


1 ' 1 ' 


D-DMRG /\ 


U = 2.2D - 


NRG 




- 


— // 




' \ _ 


-^rJ\ , 1 


, 1 V 


1 , 




^ 1.5 

3, 
q1 

Q 1 

0.5 



' 1 ' 1 


1 




ri = 0.035 D 


U = 2.6D - 


// V 

- j/ 


J 




' A 

/ \\ - 


^-J \ , 1 




y 1 , iV: 



^ 1.5 

3, 

Q 1 
0.5 




' 1 ' 1 


... 


-■ 


1 ' 1 ' 








U = 2.8D - 


// X 
// \ 

— I/ \ 






y \\ - 
f \\ 

If A ~ 


-"J 1 , 1 






1 , 1 



-2-10123 
CO/ D 

Figure 13: (Color online) Spectral densities p{i-li) in the 
metallic regime close to the metal-to-insulator transition. 
The dashed lines show results of an NRG-based DMFT 
approach. For U = 2.6D, the dark solid lines show the raw 
D-DMRG data broadened with r] = 0.035/? around the inner 
Hubbard band edges, for discussion see Sec. IIIB2 



onwards by background fits as shown in Fig. 17. The re- 
sulting values and their uncertainties are included in Fig. 
7. Their extrapolation corroborates that the pseudogap 
evolves continuously to the insulating gap at Uc2 ~ 3D 
where the metallic solution becomes unstable towards the 
insulating solution. 

The well pronounced pseudogap in the metallic solu- 
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tions close to the transition reveals clearly that the quasi- 
particle peak becomes isolated from the Hubbard bands 
at higher energies. We demonstrated in Fig. 2 in Ref. 10 
that the spectral densities of the metallic and of the co- 
existing insulating solution almost coincide at higher en- 
ergies. This means that except for the quasiparticle peak 
and the inner edges of the Hubbard bands the insulating 
and the metallic solutions arc almost identical. The two 
phases differ in their behavior at small and intermediate 
energies. This phenomenon has been baptized separa- 
tion of energy scales. The assumption of its validity is 
at the basis of the numerical projective self-consistent 
method^'*"'"''"' and of Kotliar's analytical Landau theory. 
Hence, the fact that our all-numerical treatment indeed 
displays the separation of energy scales represents a valu- 
able confirmation of this important hypothesis. 

It is very instructive to compare our findings for p{uj) 
with the NRG data from Ref. 4(J. The quantitative com- 
parison is depicted in some of the panels in Fig. 12 and 
in Fig. 13. The agreement is very good in the close 
vicinity of the quasiparticle peak, i.e., at low energies. 
But significant deviations are discernible for intermedi- 
ate and higher frequencies. The Hubbard bands obtained 
by the D-DMRG are sharper compared to the NRG re- 
sults. Hence their precursors become visible already at 
lower values of U, see Fig. 12. Furthermore, the Hub- 
bard bands in D-DMRG do not have significant tails at 
higher energies. This difference stems from the broad- 
ening proportional to the frequency which is inherent to 
the NRG algorithm. '^'""^'"'^ Similar differences were 
also reported in Ref. 32 where perturbative results are 
compared with NRG data. 
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Figure 14: (Color online) Part of the spectral density p{io) and 
the associated bath parameters y„ calculated self-consistently 
for various numbers L of fermionic sites. Left panel: U = D. 
Rigtit panel: U = 2D. 



In order to investigate a possible influence of finite-size 
effects, we consider the self-consistent spectral densities 
and the associated bath parameters 7„, cf. Eq. (14), for 
various sizes of the bath. Fig. 14 shows the results for 
L = 120, 160 and 240 fermionic sites well in the metallic 
regime ioT U = D and U = 2D. If the finite broadening 
of the D-DMRG raw data is chosen according to Eq. (19) 
the spectral densities deconvolved by LB arc not affected 
significantly by the number of bath sites. We observe nei- 
ther qualitative nor significant quantitative differences. 
This observation applies also to the self-consistently de- 
termined bath parameters 7„. Above n « 80 almost all 
bath parameters take the same value which determines 
the band width, see App. A. 

We conclude that a system size of 80 fermionic sites is 
already sufficient to capture quantitatively the features of 
the spectral properties well in the metallic regime. Nev- 
ertheless, a small improvement in the pinning criterion 
can be observed if more fermionic sites are considered, 
see Ref. 15). This improvement arises from a better ap- 
proximation of the thermodynamic limit which makes the 
dcconvolution using the LB algorithm more robust. But 
the system size cannot be chosen too large in practice 
because the D-DMRG calculation will acquire system- 
atic errors for large systems and given computational re- 
sources. Our D-DMRG calculations were all done with 
128 states, partly checked by runs with 256 states. But 
we have not noticed significant deviations between the 
runs with 128 and those with 256 states. 



2. Sharp Inner Side Peaks 

For strong interactions {U > 2D) close to the metal- 
to-insulator transition the metallic solutions show sharp 
peaks at the inner edges of the Hubbard bands, see 
Fig. 12 for U ^ 2D and Fig. 13. These peaks be- 
come sharper and sharper on increasing the interaction 
U. Note, that their height cannot (and does not) exceed 
the maximum value of the non-interacting DOS because 
the interacting DOS is determined from the Dyson equa- 
tion (4). It implies that p{lu) results from p"(ijj) by a shift 
in frequency due to *HeE and by a Lorentzian broadening 
due to UmS so that p{u!) is bounded by the maximum 
of p'^{uj). The sharp features do not represent diverging 
singularities. 

In order to properly resolve the sharp side peaks as well 
as the narrow quasiparticle peak, we choose Acui = rji = 
0.Q5D and partially Aw^ = rji = 0.025Z? in the relevant 
frequency regions. Furthermore, we place the computed 
points of the raw data at the frequency of the maximum 
of the side peaks. 

In Fig. 15 we show that the sharp peaks at the inner 
band edges are not an artifact of the deconvolution al- 
gorithm since a large variety of deconvolution schemes 
essentially provides the same peak. Furthermore, the 
side peaks are already discernible in the raw data at 
small broadening (dark solid lines) shown in Fig. 13 for 
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Figure 16: (Color online) Part of the spectral density p{u}) and 
the associated bath parameters 7„ calculated self-consistently 
for various numbers L of fermionic sites. Left: U = 2.4_D. 
Right: (/ = 2.6D. 



Figure 15: (Color online) Comparison of various deconvolu- 
tion schemes for U = 2.6D. The solid line with filled circles 
depicts the raw data. MI refers to linear matrix inversion 
for which we assume that the DOS consists of a set of S- 
functions [spikes, MI(S): circles] or that it is piecewise con- 
stant [histogram, MI(H): squares] or that it is piecewise linear 
and continuous [MI(C): triangles], for details see Ref. 23. LB 
stands for the continuous positive result obtained by the non- 
linear LB deconvolution. Upper panel: complete frequency 
interval of the symmetric DOS. Lower panel: details at lower 
frequencies. 



U = 2.6D. 

We checked that different initial hybridization func- 
tions led to the same self-consistent solutions after suffi- 
cient iterations so that we exclude that the side peaks 
are artifacts of the self-consistency procedure. More- 
over, Nishimoto et al. found very similar side peaks by a 
completely independent implementation of DMFT by D- 
DMRG."^ Hence we conclude that the sharp side peaks 
resolved in our D-DMRG approach are part and parcel 
of the metallic spectral density of the half-filled Hubbard 
model close to the metal-to-insulator transition. 

We repeated the investigation of the possible influ- 
ence of the system size in Fig. 14 close to the metal- 
to-insulator transition. In Fig. 16 we compare spectral 
densities for U = 2AD and U = 2.6D and their bath 
parameters 7„ determined self-consistently for 120, 160, 
and 240 fermionic sites with Auji = r]i = O.ID. Smaller 
values for Aw and r] turned out not to be feasible for 
the large system sizes due to limited computational re- 
sources. Thus, we accept a slight loss in the resolution of 



the side peaks compared to the spectral densities shown 
in Fig. 13 where calculations with Acj.; = jy^ = 0.0251? 
were performed for frequencies around the side peak po- 
sitions. 

We see in Fig. 16 that the precise position and the 
shape of the side peaks depend slightly on the size of the 
bath considered in the numerical calculations. The other, 
less sharp features, like the central peak and the Hubbard 
band are very robust against the changes of the bath 
size. It is very instructive to look at the bath paramters 
7„. For U — 2 AD they do not change on passing from 
L = 120 over L = 160 to L = 240. The oscillations 
visible from smaller values of ri ^ 80 are essentially the 
same for the three system sizes. They do not extend to 
larger n for larger L. For U = 2.QD, however, the bath 
parameter change for larger systems. The oscillations 
extend to higher n for L = 240 than for lower L. For 
L — 240 they appear to have converged for large n so that 
we assume that increasing L further would only extend 
the region of almost constant 7„. 

We conclude, that sufficiently large baths are essential 
for a quantitative description of the spectral densities and 
of the sharp features at the inner edges of the Hubbard 
bands in particular. This is especially important close 
to the mctal-to-insulator transition. Our results confirm 
the existence of the sharp side peaks and provide their 
accurate determination for all U < 2.2D. For U > 2.2D 
the possibility of small inaccuracies of the order of the 
deviations shown in the upper panels of Fig. 16 must be 
kept in mind. 

The only previous evidence for sharp side peaks in 
numerical investigations were weak shoulders in NRG'^'' 
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and quantum Monte Carlo (QMC)"^^ data. The inherent 
broadening in the NRG calculations and the necessity 
to continue the QMC data from imaginary times to real 
frequencies have both a tendency to smear out sharp fea- 
tures away from the Fermi level, cf. Ref. !). Hence we 
tend to interprete the shoulders in Refs. 4() and M as 
broadened remnants of the sharp peaks resolved in the 
present D-DMRG calculation. 

Besides the numerical approaches semi-analytical ap- 
proaches, the iterated perturbation theory (IPT) '^'"' and 
the non-crossing approximation (NCA), display simi- 
lar behavior of the spectral densities near the inner Hub- 
bard band edges. The IPT results show a substantial 
accumulation of spectral weight at the inner band edges 
which differs very much from our findings: the weight 
is much larger and the features are by far not as sharp 
as those presented here. They carry even more spectral 
weight than the quasiparticle peak and they are accom- 
panied by spike structures over the entire Hubbard bands 
in contrast again to our findings. In view of these differ- 
ences and because IPT is an approximate diagrammatic 
approach we conclude that the IPT features do not cor- 
respond to the ones presented here. 

The NCA data in Ref. 52 is qualitatively much closer 
to our data as far as the side peaks are concerned. Since 
the NCA violates the Fermi liquid properties for temper- 
atures lower than the Kondo scale ' * reliable conclusions 
can only be drawn for sufficiently large temperatures. 
Moreover, the peaks are found in the insulating phase so 
that it is an open question whether the NCA findings are 
related to the D-DMRG side peaks at zero temperature 
in the metallic regime. 

In order to pave the way to a physical interpretation 
of the side peaks we determine its weight S. Motivated 
by the results for insulating p{uj) in Sec. Ill A 1 we fit an 
approximate square-root onset multiplied by a quadratic 
polynomial to the Hubbard band. The excess weight is 
identified with the spectral weight S of the sharp side 
peak. This procedure is illustrated in Fig. 17. 
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Figure 17; (Color online) Determination of the side peak 
weight S illustrated for U = 2.5D. 
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Figure 18: Weight S of the peaks at inner Hubbard band 
edges as a function of the quasiparticle weight Z. 

Inspection of Fig. 13 suggests that S vanishes for 
U — > Uc2- This conclusion is supported by the fact that 
the insulating solutions do not display a similar feature, 
see Sec. Ill A 1. Hence we presume that the side peak in- 
volves quasiparticles from the central peak. The weight 
in this central peak is given by Z which vanishes linearly 
for U Uc2, see Sec. IIIB4. So we plot the results for 
S as function of Z in Fig. 18 including a point S" = for 
Z = 0. The error bars result from the numerical difficulty 
to resolve this sharp feature and from the analytical diffi- 
culty to separate it from the background of the Hubbard 
bands. 

In view of the errors the interpretation of Fig. 18 is dar- 
ing. But our guiding hypothesis is that the surprisingly 
sharp side peak is the signature of a composite excitation 
involving several fermionic excitations. ' " A single quasi- 
particle alone cannot be at the origin of a sharp feature 
because it has a finite dispersion and a finite life-time so 
that it engenders only broad features in local quantities 
away from the Fermi level. 

One possible assumption is that the composite excita- 
tion comprises three quasiparticles (Note, that the total 
statistics must be fermionic). But they need to be bound 
or anti-bound to account for the very sharp structure 
which tells us that the composite entity must be almost 
immobile. Otherwise it would give rise only to a broad 
structure in local spectral densities, which correspond to 
sums over the entire Brillouin zone. A plausible second 
scenario is that a single quasiparticle is combined with 
a collective mode, for instance a particle-hole pair as it 
occurs when a spin is flipped. 

With these ideas in mind, we want to use Fig. 18 only 
to decide whether one, two or three quasiparticles are in- 
volved. The quasiparticle weight Z measures to which ex- 
tent a local fermionic creation operator generates a true 
quasiparticle including its polarization dressing. So we 
presume that the contribution of a single quasiparticle 
gives rise to a linear behavior S (x Z while two quasi- 
particles lead to S oc Z'^ , and three quasiparticles imply 
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S <x. . Hence we find in the linear fit in Fig. 18 a liint 
that only one quasiparticle contributes to the side peaks. 

Additionally, we presume that an immobile collective 
mode is involved. Note, that due to the scaling (2) in 
DMFT collective modes are as such immobile because 
their motion requires at least two single-particle propa- 
gators. The nature, however, of the collective mode in- 
ducing the composite feature is still an open issue which 
we address in ongoing investigations. 



3. Self-Energy 

In the metallic phase we use again Eq. (26) to ob- 
tain the self-energy. The results are plotted in Fig. 19. 
The real and imaginary part of the self-energy display 
a typical low-frequency behavior of a local Fermi liq- 
uid: The real part shows a linear behavior with a nega- 
tive slope which steepens as U increases. In this region, 
the imaginary part of shows a quadratic behavior, 
3mE(ti') cx —uj^, as expected.- The quadratic behavior 
of 3mS is not perfectly retrieved, see inset of Fig. 19. 
For low values of U it seems that 3mE is proportional 
to a higher power than for w — 0. For large inter- 
action values, for instance U = 2.6D, the computation 
of the self-energy via (26) even leads to negative spectral 
densities (positive nmE for the retarded self-energy) vio- 
lating causality. This is a clear indication for the limits 
of the accuracy of the present calculation. But we stress 
that these difficulties are minor and do not affect the self- 
consistency because they occur only when the self-energy 
is computed from the full propagator G. 

The two-peak structure in the imaginary part of S(w) 
is closely related to the three-peak structure of the spec- 
tral density given by lower Hubbard band, central Kondo 
peak, and the upper Hubbard band. According to 
Eq. (26) the self-energy has a peak if the real and the 
imaginary part of G{uj) are small. For the metallic solu- 
tions, this is the case in the frequency regime between the 
Hubbard bands and the quasiparticle peak, i. e. within 
the pseudogap. 

This relation has been previously analyzed.""^''''"'''' ' It 
turns out that the strong peaks in the self-energy are lo- 
cated at position oc ±^/ZD, where Z is the quasipar- 
ticle weight, see next subsection. Our data agrees very 
well with the analytical arguments; the factor of propor- 
tionality is about \/2. We stress that ^ on ap- 
proaching the metal-to-insulator transition whereas the 
sharp inner side peaks remain at the inner band edges. 
Hence the two peaks in the imaginary part of the self- 
energy are not directly related to the sharp side peaks. 
It should be pointed out that the two-peak structure of 
the self-energy contradicts the findings in Ref. •')5. We 
presume that the analytical argument misses some effect 
which requires further research. 




Figure 19: (Color online) Real and imaginary part of the self- 
energy E(tj) in the metallic regime. Inset: 3m S(a;) on a 
larger scale. 



4- Quasiparticle Weight 

The quasiparticle weight Z is the characteristic quan- 
tity of the metallic phase. Crudely speaking Z quantifies 
how much of the spectral weight is found in the central 
quasiparticle peak, cf. Fig. 13 and 12. In other words, Z 
is the matrix element between the fermionic creation op- 
erator applied to the ground-state and the dressed quasi- 
particle excitation including its polarization cloud. For 
local self- energies- it is given by 



z-^ = \- a^s(w)L 



(31) 



Because the self-energy is not directly accessible in our 
approach we prefer to calculate Z using the derivative of 
the Dyson equation (4) implying 



(32) 



where d^G^{uj)\^=Q = "^/d^ has been used. By Eq. (32) 
we avoid the numerical errors occurring if the self-energy 
is determined via (26), see Sec. IIIB3. 

The derivative dLjG{uj)\u,=o in Eq- (32) is real because 
the imaginary part displays an extremum at a; = 0. The 
derivative of *He G{uj) is obtained reliably from the coef- 
ficient a of the fit function 



/meG(w) = aw + buj^ 



(33) 
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Figure 20: (Color online) Quasiparticle weight Z as a, function 
of U compared to perturbative results '" (PT) up to 0(11''') 
and to NRG data.^" 

where the absence of even orders due to particle-hole 
symmetry has been taken into account. The fit is done in 
the intervals around the origin which decrease for U 
Uc2 roughly linearly; for U = 2D it is [-0.1D,0.1D]. 
The derivative cannot be determined robustly as ratio 
of differences because 9{eG{uj) displays weak oscillations 
resulting via the Kramcrs-Kronig transformation from 
weak wiggles in 3mG(w) due to the deconvolution proce- 
dure, see for instance the lowest panel in Fig. 13 around 

LU = D. 

In Fig. 20 the quasiparticle weight Z is depicted as a 
function of U. Starting from Z ^ \ a±U — 0, Z decays as 
U increases. Our results agree excellently with perturba- 
tion theory where it is applicable U ^ D. The NRG data 
is lower for all values of U. We deduce that the broaden- 
ing of the spectral weights proportional to w, inherent to 
NRG, has a non-negligible influence on the quasiparticle 
weight. Qualitatively, however, the behavior of Z{U) is 
very similar. 

From the trend of Z{U) one expects that Z vanishes 
continuously at a finite critical value Uc2- We cannot 
approach this value arbitrarily closely because the reli- 
able resolution of the sharper and sharper quasiparticle 
peak would require a better and better resolution. This 
in turn cannot be achieved with a reasonable amount of 
computing resources because Eq. 19 must be respected. 
For this reason, we have to resort to an extrapolation to 
determine Uc2 quantitatively. 

The behavior of Z{U) close to Uc2 was already in- 
vestigated previously analytically by the projected self- 
consistent method (PSCM)^''* and by the linearized 
DMFT. ■'' Both methods assume the separation of en- 
ergy scales and predict a linear behavior of Z close to 
Uc2- Since our all- numerical spectral densities strikingly 
confirm the separation of the energy scales^" we deter- 
mine Uc2 by fitting 

/z(f/) = ai(C/-C/c2)+a2(f/-C/c2)' (34) 
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to the data points for U G [2, 2.8]!?. We obtain a crit- 
ical value of Uc2 = (3.07 ± 0.1)13 which agrees excel- 
lently within the error range with most of the previous 
results, ''-'' *''''-^^''-"^'**^'***''^^ including the hnearized DMFT"'" 
and the PSCM."'^ For comparison, we include an overview 
of these critical values in Tab. L 



(7ci/D 




method 


2.38 ±0.02 


3.07 ±0.1 


D-DMRG (CV, direct inversion) 


2.225 ± 0.025 


3.05 ±0.05 


D-DMRG (CV, variational)-''-'^ 


2.39 ±0.02 


3.0 ±0.2 


D-DMRG (Lanczos) " 


2.39 




PT-QMC'"' 


2.39 


2.94 


NRG-'"'^' 




3 


hnearized DMPT"'" 




2.92 ±0.05 


PSCM'^ 



Table I: Comparison of the values for Uci and Uc2 found by 
various approaches. 



5. Momentum-Resolved Spectral Densities 

For the metallic phase, we present spectral densities 
p{uj, efe) as defined in Eqs. (28,29) which arc resolved ac- 
cording to their bare single-particle energy . In the con- 
text of the DMFT, this procedure provides information 
on the momentum dependence of the spectral proper- 
tics. Figs. 21 and 22 display how the momentum-resolved 
spectral densities evolve on approaching the metal-to- 
insulator transition. 

Three main features can be discerned. First, the dom- 
inant diagonal line of highest intensity with Whigh oc e/c 
catches our attention in all six panels. This is where the 
Fermi liquid picture manifests itself. " Without interac- 
tion a quasiparticle can be excited at o^high ~ Cfc- In the 
presence of interaction there are still quasiparticlcs which 
behave like fermions above the Fermi level or holes below 
it. But these excitations are dressed by clouds of particle- 
hole pairs. This dressing has two effects. One is that the 
bare fcrmionic operator has only an overlap reduced by 
Z with the creation of the quasiparticle. The other effect 
is that the dressed quasiparticlcs are heavier than their 
bare counterparts. This means that their dispersion is 
reduced. For a momentum-independent self-energy the 
reduction of the dispersion is given again by the factor 
Z < 1. This can be seen in the impressive steepening of 
the dominant diagonal line given by Whigh = Zek upon 
increasing interaction. For interaction values close to 
the instability of the metallic solution, sec lowest panel 
in Fig. 22, the diagonal line has become so steep that 
it can almost be taken as a vertical one. This implies 
that the quasiparticle has become almost local and non- 
dispersive. This is the mechanism leading to the so-called 
heavy fermions. 

We point out that the white stripes at w = in most 
of the panels in Figs. 21 and 22 and the white stripes at 
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Figure 21: (Color online) Spectral densities p{u,ek) in the 
metallic regime for U = 0.5-D (top), U = D (middle), and 
U = 2D (bottom). 



Figure 22: (Color online) Spectral densities p{uj,ek) in the 
metallic regime for U = 2.4_D (top), U — 2.6D (middle), and 
U = 2.SD (bottom). 



small uj in the first panel of Fig. 21 result from the numer- 
ical violations of causality as discussed before, see Sects. 
Ill A 3 and IIIB3. For values of uj where this problem 
occurs the spectral densities are set to zero. 

The second main feature in the momentum-resolved 
spectral densities is the development of the precursors of 
the lower and the upper Hubbard band. At C/ = Q.5D 
no such precursors are visible, except for some broaden- 
ing of the diagonal line away from the Fermi level. But 
in the panel for U — 1.5D the precursors are already 
clearly visible. They are qualitatively very similar to the 
broad densities found in the insulating spectral densities 
in Figs. 9 and 10. This similarity becomes more and 
more pronounced as U approaches 11^2 where the metal 
changes continuously to the insulator. 

The third main feature is the occurrence of a vertical 
line of high density at the inner edges of the precursors 
of the Hubbard bands. It is barely visible for U = 2.0D, 



see last panel of Fig. 21. But it is clearly discernible 
for the larger values of the interaction displayed in Fig. 
22. This vertical line represents the momentum-resolved 
feature of the sharp side peaks at the inner band edges 
discussed in Sec. TUB 2. We re-emphasize that there is 
no discernible momentum dependence indicating a local, 
non-dispcrsivc, immobile excitation. 

We draw the reader's attention to the fact that the 
vertical lines display an interesting dependence of their 
weight on the momentum. The high density is found for 

> for efe < and vice versa {uj < for Ck > 0). 
This implies that the immobile composite entity is ex- 
cited upon adding a bare fermion {uj > 0), but for mo- 
menta corresponding to occupied states (e^ < 0) and vice 
versa. We deduce that the addition of the bare fermion 
must induce the creation of a quasiparticle hole, cf. Sec. 
IIIB2. Further work to elucidate this puzzling behavior 
is called for. 
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IV. SUMMARY 

The present work constitutes the comprehensive anal- 
ysis of the transition from a metal to an insulator driven 
by strong repulsive interaction. The framework of this 
investigation is the dynamic mean-field theory (DMFT). 
We evaluate the necessary self-consistency equations at 
zero temperature with a numerical impurity solver based 
on dynamic density-matrix renormalization (D-DMRG). 
The D-DMRG computes a frequency-dependent correc- 
tion vector which ensures a high reliability of the re- 
sults at all energies, not only close to the Fermi level. 
Thereby, we are able to provide spectral densities in the 
metallic and in the insulating phase with unprecedented 
resolution. Such densities are relevant to photoelectron 
spectroscopy of transition metal compounds and other 
narrow-band systems. 

The calculations are done for a semi-elliptic non- 
interacting density of states (DOS). All the important 
technical details of the intricate combination of D-DMRG 
and DMFT arc given in detail so that our findings may 
be reproduced independently. 

The local, momentum-integrated, spectral density in 
the insulating phase is governed by the lower and the 
upper Hubbard band. Both are featureless and display 
square-root singularities at the inner band edges. The 
insulating gap A vanishes continuously for J7 — > Ud = 
(2.38 ± 0.02)£'. A slight curvature indicates that A cx 
{U - C/ci)« with C = 0.72 ± 0.05. 

For U < Uc2 = (3.07 ± 0.1)1) the insulating solution 
is only metastable. The phase with lower ground-state 
energy is the metallic one.^" The corresponding spectral 
densities show three important features for significant in- 
teraction values: (i) the central quasiparticle peak, (ii) 
the precursors of the lower and the upper Hubbard band, 
and (iii) the sharp side peaks at the inner band edges of 
the Hubbard bands. 

The DOS of the central peak at the Fermi level is 
pinned to its non-interacting value due to Luttinger's 
theorem. The central peak is built from heavy fermions 
or holes which disperse only weakly due to their strong 
dressing. They are the characteristic signature of a Fermi 
liquid. 

The Hubbard bands signal that the metal evolves to- 
wards an insulator on increasing interaction. Our find- 
ings, both the earlier ones'" and the present ones, clearly 
show that the metallic spectral densities evolve continu- 
ously to the insulating one upon U Uc2- This continu- 
ity holds in the sense of an integral norm, i.e., considering 
the distribution of spectral weight, not in the sense of an 
absolute value norm which jumps. For instance, the DOS 
at the Fermi level is 2/{ttD) in the metallic regime and 
zero in the insulating one. But the corresponding weight 
vanishes linearly for U Uc2- The pseudogap devel- 
ops simultaneously to the insulating gap. These findings 
numerically confirm the assumption of the separation of 
energy scales.^*' '"" 

The sharp side peak at the inner edges of the Hub- 



bard bands in the metallic regime are the signature of 
a complex composite excitation which is not yet under- 
stood. We found that its weight scales with the quasi- 
particle weight which supports our hypothesis that the 
side peaks stem from an entity made from a quasiparti- 
cle and some collective, immobile mode. The analysis of 
the momentum-resolved spectral densities shows in ad- 
dition that the involved quasiparticle is a hole if a bare 
electron is added. Conversely, it is a particle if a bare 
hole is created, i.e., an electron is taken out. 

Ongoing research is devoted to the nature of the com- 
posite excitation which generates the sharp side peaks. 
To this end, the response functions, which correspond to 
various collective modes, need to be investigated. 

Other interesting issues concern the qualitative and 
quantitative dependence of the important features found 
on the non- interacting density of states. The changes 
introduced by doping certainly will be of particular in- 
terest. A first study, based on Lanczos-DMRG, of the 
effects of doping has been published very recently. ""^ 
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Appendix A: CALCULATION OF THE 
CONTINUED FRACTION COEFFICIENTS 

Once a propagator g{Lu) is known on a dense mesh 
of frequencies (about 10^ points), one is able to deter- 
mine its continued fraction representation to an arbitrary 
depth 

9{^) = '—yl (Al) 

c-ao '-^^ 

bj — ai 

w - a2 - • • • 

very precisely. Numerically, only integrations are in- 
volved. 

For the iteration, we define recursively the resolvents 
Po(t^) := g(t^) , (A2a) 
p„+i(a>) = w - a„ J—. (A2b) 

Obviously, the continued fraction representation of p„(ti;) 
starts with a„ and &„. It is shifted by n relative to the 
continued fraction of _g. An expansion in 1/cj reads 

p„(c.) = ^(l + ^-KO(c.-2)). (A3) 
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Comparing this expression with the expansion in 1 jio of 
the Hilbert representation of p„(a;)""'^^ leads to 

oo 

bl^-'^ J duj3mpn{uj) , (A4a) 

— OO 

an = --^ dujuj3mpn{uj)- (A4b) 
T^bl J 



The iteration of the above equations for n ~ 0,1,2,... 
yields the wanted coefficients. The iterations are stopped 
once all the coefficients for the numerically tractable fi- 
nite chain of L fermionic sites are determined, i.e., for 
n ~ L — 3. 

In the case of a particle-hole symmetric model, the 
spectral density is symmetric. Consequently, all odd mo- 
ments and thus all the coefficients a„ (n = 0, 1, . . .) van- 
ish in the above recursion scheme. So for the model under 
study we need only to consider the equations for the 6„. 

The asymptotic behavior of the continued fraction co- 
efficients for n — > cxD is determined by the band width 
of the spectral density and by the singularities at 
the band edges.'"' We exploit the relation for the band 
width 2B, if the support of the spectral density is com- 
pact, which reads 



4|6„ol = 2B 



(A5) 



where 6oo := limi^oo^i- For a gap A in the middle of 
the spectral density, i.e., for the insulator, the bi are os- 
cillating for large i like 



feoc±(-l)*A/2 



(A6) 



where &oo is again given by a quarter of the band width 
which is the distance between the largest and the smallest 
band edge including the gap. 

In a strict sense we expect that the support of the spec- 
tral density of an interacting system is infinitely large. 
But it turns out that significant spectral weight is found 
only in a finite interval. The contributions outside this 
interval are negligible which has been observed also in 
other calculations."'^ 

We use the relations (A5,A6) to determine the inter- 
val in which we have to compute raw data {gi}. For 
this purpose we compute boo in each iteration of the self- 
consistency cycle and observe its development on itera- 
tion. If the value of b^o does not change significantly for 
several iteration steps, we fix the border B of the inter- 
val in which we perform the calculation of the raw data. 
This procedure allows us to avoid the time-consuming 
calculation of unnecessary data points. 

In the presence of an insulating gap A at the Fermi 
energy the above recursion scheme (A4) to determine 
the continued fraction coefficients has to be modified for 
practical reasons. Due to the antisymmetry of its real 
part the propagator in the particle-hole symmetric case 
vanishes at a; = completely pq{u! = 0) = with no 



imaginary part. So there is a pole in pi{u!) at w = 
according to (A2b) implying a finite contribution from 
a ^-distribution at a; = in 3m pi. We account for it 
separately 



oo 

5? = / doj Jmpiiuj) 

TT J 



1 



(A7) 



[-A,A] 



The recursion relation (A2b) implies that the same sce- 
nario occurs for every odd n. So we use (A7) for any odd 
n once the subscript i is replaced by „ and the subscript 
by n-i- 



Appendix B: REMOVAL OF DECONVOLUTION 
ARTIFACTS 



1. Outer Band Edges 

The raw data {gi} from D-DMRG covers only a certain 
finite interval. Hence the underlying unbroadened spec- 
tral density p{uj) can also be determined only on a finite 
interval [—B, B] {B « 3D) which is very close to the 
interval in which the raw data is available. The decon- 
volution leads to an abrupt increase of the deconvolved 
p{llj) at the assumed band edges at B, see Fig. 23. This 
increase is definitely an artifact. In general, we expect 
that the spectral density p{uj) continues to decrease for 
I a; I — > oo. This decrease is a gradual one, i.e., there is no 
band edge in the strict sense for the correlated system. 
But all results obtained so far by us and others"'' point 
into the direction that the contributions beyond a certain 
interval can be safely neglected. 




Figure 23: Lines: Generic deconvolved spectral density close 
to the band edge at B. The dashed line marks the artifact 
which is removed by setting p{uj > tJmin) ~ p(i^min). The left 
band edge is treated accordingly. 

The artificial increase is removed by replacing p{uj) for 
frequencies beyond cUmin where the first minimum occurs 
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by the minimal value p{iL>min)- This is indicated by the 
dashed line in Fig. 23 which is replaced by the solid curve. 

2. Inner Band Edges of the Insulator 

The situation close to inner band edges stemming from 
a gap is similar to the one close to outer band edges. 
The LB dcconvolution of the D-DMRG raw data in the 
region, where the insulating solution displays a gap. leads 
to small artifacts as shown in Fig. 24 by the dashed lines. 




Figure 24: Lines: Generic deconvolved spectral density in the 
region of the gap A of the insulating phase. The dashed line 
marks the artifacts which are cut out by setting p(uj) to zero 
for io G [-A5, Aa]. 

The first remedy, which one may consider, is to deter- 



mine the gap A''*'*' of the SIAM by conventional static 
DMRG and to use this value as the gap of the lattice 
problem in DMFT. But we found out that even an ex- 
trapolation of this static gap A'^''** to the thermodynamic 
limit L — > 00, as proposed in Ref. 27, cannot be used in 
our iterations. No stable self-consistent solution in the 
sense of the criteria given in Sec. II C 4 could be obtained. 

The successful treatment of the artifacts consisted in 
removing only the obviously artificial features which are 
depicted by dashed lines in Fig. 24. The solid lines are 
kept because they represent real spectral weight. This is 
done by starting from the value of the static gap A'^'''* 
adjusted by a correction parameter S > 

p{u}) = Vcj e [-A'^'^' + iJ, A'^'^'-.J]. (Bl) 

The value 6 is chosen carefully, so that on the one hand 
all dcconvolution artifacts are removed and on the other 
hand the proper spectral weight of the Hubbard bands 
stays untouched. In order to determine the appropri- 
ate value of 5 it is helpful to resolve the inner Hubbard 
band-edges in detail, which can be achieved by small val- 
ues of r]i in this region. Note, that the adjusted interval 
[— A'^'^* _l_ ^stat _ j-jg smaller than the interval 

[—A, A] defined by the proper gap A as long as all dccon- 
volution artifacts are removed. The artifacts would spoil 
the DMFT self-consistency iteration because the contin- 
ued fraction coefficients could not be evaluated properly. 

It turns out that 6 does not need to be changed in 
each iteration of the DMFT self-consistency cycle. It is 
sufficient to set it once for each value of U. 
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